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ABSTRACT 

Redshifted 21cm signal is a promising tool to investigate the state of intergalactic 
medium (IGM) in the Cosmic Dawn (CD ) and Epoch of Reionization(EoR). In our 
previous work (|Shimabukuro et al. 20l'd) . we studied the variance and skewness of 
the 21cm fluctuations to give a clear interpretation of the 21cm power spectrum and 
found that skewness is a good indicator of the epoch when X-ray heating becomes 
effective. Thus, the non-Caussian feature of the spatial distribution of the 21cm signal 
is expected to be useful to investigate the astrophysical effects in the CD and EoR. 
In this paper, in order to investigate such a non-Caussian feature in more detail, we 
focus on the bispectrum of the 21cm signal. It is expected that the 21cm brightness 
temperature bispectrum is produced by non-gaussianity due to the various astrophys¬ 
ical effects such as the Wouthysen-Field (WE) effect. X-ray heating and reionization. 
We study the various properties of 21cm bispectrum such as scale dependence, shape 
dependence and redshift evolution. And also we study the contribution from each 
component of 21cm bispectrum. We find that the contribution from each component 
has characteristic scale-dependent feature. In particular, we find that the bulk of the 
21cm bispectrum at 2 ; =20 comes from the matter fluctuations, while in other epochs 
it is mainly determined by the spin and/or neutral fraction fluctuations and it is ex¬ 
pected that we could obtain more detailed information on the ICM in the CD and EoR 
by using the 21cm bispectrum in the future experiments, combined with the power 
spectrum and skewness. 

Key words: cosmology: theory — intergalactic medium — Epoch of Reionization 
— 21cm line 
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1 INTRODUCTION 

After the recombination of protons and electrons in the pri¬ 
mordial plasma, the “dark age” begins and continues un¬ 
til the first object is formed through the gravitational in¬ 
stability in dark matter halos. The universe dawns by the 
lights from the first stars and astrophysical processes start 
to play important roles in the intergalactic medium (IGM) . 
( Yoshida et al. 20061 : iFialkov et al. 201^ : IVisbal et al. 20 f2 : 
ISantos et al. 2 OIII I. Among others, Wouthuysen Field (WF) 
effect due to the Lyman alpha (Ly-a) radiation emitted 
from the primordial stars is expected to be the first impor¬ 
tant process for IGM, which couples the spin temperature 
of neutral hydrogen to the color temperature of Ly-a ra¬ 


diation (I Wouthuysen 195^ . In the typical environments of 
IGM, color temperature of Ly-a radiation field couples with 
kinetic temperature of IGM because of the large Ly-a scat¬ 
tering rate. Therefore, the spin temperature also couples to 
kinetic temperature. The second astrophysical process which 
is related to the the rmal history of IGM after the WF ef¬ 
fect is X-ray heating llPritchard fc Furlanetto 2007l l. The X- 
rays are thought to be emitted from su ch sources as galaxies 
and X-ray binaries llFialkov et al 20141 ) , and kinetic temper¬ 
ature of IGM increa ses dramatically jMesinger et al 20121 : 
IChristiari_et_ahj01^ WFinallv. the “Epoch of Reionization 
(EoR)” ( Fan et al. 2006l l follows these processes when the 
gas density of IGM decreases enough so that recombination 
process becomes inefficient. This is the epoch when neutral 
hydrogens star t to be ionized by UV r adiation photons from 
early galaxies dLoeb &: Barkana 200ll L 

The redshifted 21cm line from neutral hydrogen due 
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to the hyperfine transition is suitable for studying ther¬ 
mal and ionized states of IGM as well as the first ob¬ 
jects in the dark age and EoR. dFur lanetto et al 200^ : 
IPritchard fc Loeb 20ld : IShimabukuro et al.2014l h One of 
the statistical methods to subtract the information about 
the physical state of IGM at those epochs from 21cm 
line is the power spectrum a nalysis of brightness temper¬ 
ature (iFurlanetto et al 20061: Pritcliarf_&_^F\irlanetto_2007 ; 


Santos et al. 20081: iBaek et al. 20101 : iMesinger et al. 2013 


Pober et al b. 2014h . On-going radio interferometers, such 


as Low Frequency Array (LOFAR) dRottgering et al. 2013ll . 
Murchison Wide Field Array (MWA) ( Tingav et al. 2012 1 
and Probing the Epoch of Reionization (PAPER) 
dPober et al a. 2014l L have started observation and serve as 
a “prototype” of future high-sensitivity experiments. Al¬ 
though the sensitivities of the on-going experiments are in¬ 
efficient to form images of the distribution of neutral hydro¬ 
gen using the redshifted 21cm line from the cosmic dawn 
and EoR, th e power spectrum o f the 21cm signal would 
be detected dMesinger et al. 2013l L In particular, it is ex¬ 
pected that the 21cm power spectrum for each redshift 
(z = 10 ~ 30) can be measured by the Squa re Kilometre 

Array (SKA) with high accura£^_JCajilli_201^ _ 

In our previous work (IShimabukuro et al. 201^ . we 
gave an interpretation to the time evolution of the 21cm 
power spectrum and we find that the size of skewness is 
sensitive to the epoch when X-ray heating becomes effec¬ 
tive. Other work also reports the impact of spin temperature 
fluctuations on the skewness dWatkinson fc Pritchard 2015l ~l 
and there is a work which focuses on the redshift distor¬ 
tion as the indicator of the epoch where X-ray heating is 
effective dFialkov et al. 2015l ~l. Herein we extend these pre¬ 
vious works by considering the the bispectrum to inves¬ 
tigate the dependence of the skewness on scales because 
skewness is an integral of the bispectrum with respect to 
the wave number (s ee Appendix A). In a different work 
dYoshiura et al. 2015f ) we have already estimated errors from 
the thermal noise of detectors in estimating the bispectrum 
and we found that the 21cm bispectrum would be detectable 
at large scales at fc ^ O.lMpc”^ even by the current detec¬ 
tors on, such as, the MWA and PAPER. Furthermore, the 
21cm bispectrum would be detectable even at small scales 
with the SKA dYoshiura et al. 2015l L and therefore the 
study of the bispectrum is timely and well mo tivated. Some 
previous works stud i ed the 21cm bispec t rum dCoorav 20051 : 
Pillepich et al. 20061 : iGoorav et al. 20081 : iMunoz et al. 20151 : 
2). These works, however, mainly focus on the bispectrum 
as a measurement of primordial non-gaussianity in mat¬ 
ter fluctuations. In our study, we instead focus on non- 
Gaussianity in 21cm fluctuations induced by astrophysical 
effects, whose size is expected to be larger than that in 
matter fluctuations. (For other probes o f non-Gaussianity 
such ans Minkowski function als, see dGleser et al.200a : 
iLee et al. 20081 : iHong et al. 20141 : 13). 1 


2 FORMULATION AND SET UP 

2.1 Formulation for the 21cm bispectrum 

A fundamental quantity of 21cm line is the brightness 
temperature, which is described as the spin tempera- 


ture offsetting from G MB temperature, given by (see, e.g, 
dFurlanetto et al 2006l H 


6Tb{u) = 


-(1 - e“^“o) 


1 + z 
27XH{1 + Sm) 

1 + z 0.15 
10 
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dVr/dr + H 

/rtbh^ 

V 0.023 


T 

1 -^ 

Ts 


[mK]. 


( 1 ) 


Here, Ts and Tj respectively represent gas spin temperature 
and CMB temperature, is the optical depth at the 21cm 
rest frame frequency uq = 1420.4 MHz, xa is neutral fraction 
of the hydrogen gas, 5m(x, z) = p/p— 1 is the evolved matter 
overdensity, H{z) is the Hubble parameter and dVr/dr is the 
comoving gradient of the gas velocity along the ling of sight. 
All quantities are evaluated at redshift z = uq/u — 1. 

Let us focus on the spatial distribution of the bright¬ 
ness temperature. The spatial fluctuation of the brightness 
temperature can be defined as 


52i{:>i) = 5n{:x.) - {Sn) ( 2 ) 

where {STb) is the mean brightness temperature obtained 
from brightness temperature map and (...) expresses the 
ensemble average. From this definition, we have the power 
spectrum of ^21 defined as 

(52i(k)52i(k')) = (27r)®5(k 4- k')P2i(k), (3) 

If the statistics of the brightness temperature fluctuations 
is pure Gaussian, the statistical information of the bright¬ 
ness temperature should be completely characterized by the 
power spectrum, and in the above expression for the bright¬ 
ness temperature given by Eq. ©, if both of the spin tem¬ 
perature and the neutral fraction are completely homoge¬ 
neous, the statistics of the brightness temperature fluctua¬ 
tions completely follows that of the density fluctuations Sm ■ 
However, in the era of GD and EOR, it is expected that the 
spin temperature and the neutral fraction should be spa¬ 
tially inhomogeneous and the statistics of the spatial fluctu¬ 
ations of those quantities would be highly non-Gaussian due 
to the various astrophysical effects. Accordingly, the statis¬ 
tics of the brightness temperature fluctuations would devi¬ 
ate from the pure Gaussian and it should be important to 
investigate the non-Gaussian feature of the brightness tem¬ 
perature fluctuations. Although such a non-Gaussian fea¬ 
ture can be investigated through the skewness of the one- 
point distribution functio ns as done in our previous work 
dShimabukuro et al. 201a L the scale-dependent feature has 
been integrated out in the skewness. On the other hand, 
the higher order correlation functions in Fourier space such 
as a bispectrum and a trispectrum characterize the non- 
Gaussian features and also have the scale-dependent infor¬ 
mation. Here, in order to see the non-Gaussian feature of 
the brightness temperature fluctuations <521, we focus on the 
bispectrum of ^21 which is given by 

(<52i(ki)(52i(k2)(52i(k3)) = (27r)®(5(ki-fk2-l-k3)i3(ki, k2, ks). 

( 4 ) 

In order to characterize the shape of the bispectrum in fc- 
space, we use an isosceles ansatz which is defined as fci = 
fc 2 = fc = afcs (a ^ 1/2). For examples, in case with a 1 
the shape of the bispectrum is often called as 'Squeezed type'' 
or “local type", in case with a = 1 it is called as “ equilateral 
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type'’, and in case with a = 1/2 it is called as ’’folded type”. 
Note that we relax the configuration condition because we 
calculate the bispectrum from the grid point. We regard the 
length within the range of 10% of side of the triangle we 
desire as the that of triangle. 

2.2 Calculation of the 21cm bispectrum 

In this paper, we calculate the bispectrum of the brightness 
temperature fl uctuations (21cm bispectrum) by making us e 
of 21cmFAST llMesinger et al. 20071 : iMesinger et al. 201lll . 
This code is based on a semi-analytic model of star/galaxy 
formation and reionization, and makes maps of matter den¬ 
sity, velocity, spin temperature, ionized fraction and bright¬ 
ness temperature at the designated redshifts. 

We perform simulations in a (200Mpc)® comoving box 
with 300® grids, which corresponds to 0.66 cMpc resolution 
or ~ 12.7(14.1) arcsec at 80 (127) MHz {z = 17 (10)) and 
1.07(1.19)deg® field of view at 80 (127) MHz {z = 17 (10)), 
from z = 200 to z = 8 adopting the following parame¬ 
ter set, {(,(:x,nir,Rm{p) = (31.5,10®®/A^o, 10^ K,30 Mpc). 
Here, ^ is the ionizing efficiency, (^x is the number of X-ray 
photons emitted by source per solar mass, T^ii is the min¬ 
imum virial temperature of halos which produce ionizing 
photons, and 7?mfp is the mean free path of ionizing photons 
through the IGM. In our calculation, we also ignore, for sim¬ 
plicity, the gradient of peculiar velocity whose contribution 
to the brightness t emperature is relatively small (a few %) 
dChara et al. 2014 1. We perform 10 realizations of simula¬ 
tions with different initial condition of density fluctuations 
and obtain brightness temperature maps. Then we evaluate 
the average bispectrum as 

_ I ^ 

abs[B{k)] = -jj abs[B{k)]i 

i=l 
1 ^ 

i 

Here, N is the number of realizations and k is the absolute 
value of k. 


3 RESULT 

In this section, we summarize our result for the 21cm bis¬ 
pectrum. 

3.1 Scale-dependence of 21cm bispectrum 

First, in order to see the scale-dependence of the 21cm bis¬ 
pectrum, we focus on the equilateral shape, that is, a = 1 
case for the isosceles ansatz discussed in the previous sec¬ 
tion. We plot the equilateral type bispectrum as a function of 
wave number k with 1-a sample variance for several redshifts 
[z = 10,15, 20 and 27) in Fig[T] Here, we use the normalized 
bispectrum which is given by fc®afes[i3(fc)]. z = 10 is a typical 
redshift during EoR, and z = lb and 20 are expected to be 
a transition time from CD to EoR, while 2 : = 27 is a typical 
time during CD. As you can see, the variance is relatively 
small and the cosmic variance is not so serious for the field 
size and wavenumbers we chose. From this figure, we can hnd 



k[Mpc'^] 


Figure 1. Equilateral type bispectra as functions of wave number 
at 2 : = 27 (red), 20 (cyan), 15 (green), 10 (purple). The shaded 
region associated with each line represents l-cr sample variance 
estimated from 10 realizations. 

that, except for the case with 2 = 20, the normalized bis¬ 
pectrum is almost scale-invariant for the equilateral shape. 
On the other hand, for 2 = 20, the normalised bispectrum 
has a scale-dependence as (xD. Such difference is expected 
to depend on what component gives a dominant contribu¬ 
tion to the 21cm bispectrum. As we will discuss this issue 
later, dominant component of the 21cm bispectrum at 2=20 
is matter fluctuations whereas dominant components of the 
21cm bispectrum at other redshifts are not matter fluctu¬ 
ations, but spin temperature and neutral fraction fluctua¬ 
tions. Since the bispectrum of dark matter fluctuations has 
scale dependence (larger at smaller scales due to the nonlin¬ 
ear gravitational growth), the 21cm bispectrum traces this 
scale dependence. 

Next, we show a comparison of bispectra of equilateral 
(a = 1), folded (a = 1/2) and squeezed (a = 10) types in 
Fig. [ 2 ] Here, the bispectra are plotted as functions of fcs for 
several redshifts. From this hgure, we can see that the scales 
and shapes which mostly contribute to the skewness since 
skewness is the integral of bispectra (see Appendix A). Espe¬ 
cially, smaller scales contributes to the skewness at 2 = 20 
although the bispectrum is nearly scale invariant at other 
redshifts. We also hnd that the squeezed type bispectra at 
2=15, 27 turn upwards at fc > 0.3Mpc“®. The squeezed-type 
bispectra brings information on smaller scales and matter 
huctuations are dominant at smaller scales. Thus, contribu¬ 
tion from matter fluctuations results in these upwards at 
squeezed type of bispectra. 

3.2 Redshift evolution of 21cm bispectrum 

Next, we consider redshift evolution of 21cm bispectrum. 
Before we show redshift evolution of the 21cm bispectrum, 
as reference, we show redshift evolution of temperatures in 
our model in Fig|3] We also show ionized evolution and the 
bispectra as functions of redshift for several a in Fig. [T] the 
equilateral shape (a = 1), the folded shape {a = 1/2) and 
the squeezed shape (a = 10) with fc = 1.0 Mpc“®. For the 
equilateral and folded cases, we can see two peaks located at 
around 2 = 20 and 12. These peaks can also be seen in the 
power spectrum of the brightness temperature fluctuations. 
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Figure 2. Scale dependence of bispectra for (fci : k 2 : fcs) = (1 : 1 : 1), (1 : 1 : 2)and (10 : 10 : 1), as functions of ^ 3 . 


P 2 i(k), with k ~ 1.0 Mp c~^ (see, e.g., our previous paper 
llShimabukuro et al. 201^ '). On the other hand, in case with 
the squeezed shape, three peaks appear at around z = 23, 
17, and 12. This feature is si milar to that of the power spec¬ 
trum with k ~ 0.1 Mpc“^ (|Shimabukuro et al. 201,^ 1. For 
the squeezed type, we take the parameter a to be 10 and 
this means ks = 0.1 Mpc“^. Hence, the squeezed-type 21cm 
bispectrum is expected to be described in terms of not only 
the power spectrum with larger two wave number (fci and 
k 2 in our case) but also that with smaller one wave num¬ 
ber (fcs in our case) and also it would have the information 
about the correlation between the long and short wavelength 
modes in Fourier space or local non-linearity in real space. 
The dip at z ~ 20 for the a = 10 case results in mode cou¬ 
pling between long and short wavelengths. As we show later 
in FigO the 21cm bispectrum of equilateral type at large 
scale(fc = 0.1Mpc“^) also shows a dip at ^ ~ 20 as simi¬ 
lar as the power spectrum as function of redshift and the 
a = 10 traces this feature. Therefore, we can conclude that 
the squeezed type bispectrum has information both on large 
scale and small scale. 


We will also investigate what physics cause such a cor¬ 
relation between the long and short wavelength modes in 
the 21cm bispectrum in later subsection. 



10 15 20 25 30 35 


z 

Figure 3. The redshift evolution of temperatures. We show evolu¬ 
tion of the kinetic temperature(red), the spin temperature(green) 
and the CMB temperature(blue). 


3.3 Decomposition of 21cm bispectrum 

As we saw in Eq. the fluctuations in the brightness tem¬ 
perature are contributed not only from the matte density 
field, but also from the fluctuations of the spin temperature 
and neutral fraction, aside from the gradient of peculiar ve- 
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Figure 4. {Top) ionization history in our model. {Bottom) Com¬ 
parison of bispectra of typical triangle configurations. We fix 
k = 1.0 Mpc“^ and take a = 1 (equilateral: red solid line), 
OL = ll2 (folded: green dashed line) and a = 10 (squeezed: blue 
dotted line) for the isosceles ansatz. 


locity which we neglect here. In this section, we decompose 
the bispectrum into the contributions from these compo¬ 
nents. 

We can rewrite Eq. o as, 


(5Ti,(x) = + (52 ;h(x))( 1 + (5m(x))(l + 5,,(x)), (6) 

where STb is the average brightness temperature and evalu¬ 
ated as, 


5Tb = 2TxWrj 


1 + z 
10 


1/2 


0.15 ffibC\ 
^mh?) V 0-023 y’ 


(7) 


with ITa being the volume average of xh- Here we charac¬ 
terize the contribution of the spin temperature T s by a new 
variable r; — l — T^/Ts (IShimabukuro et al. 201^ . By using 
this parameter, we can take into account a non-linear rela¬ 
tion between the spin and brightness temperatures linearly. 
The volume average of rj Is represented by rj. Note that when 
Sts ^ 1) we have. 


Sr, 


E/fs , 


( 8 ) 


Using Eq. , we can decompose the brightness temper¬ 
ature bispectrum into auto- and cross-correlation of Sm, Sx^ 
and Sr,'. 

BsTb = 

-I-(cross correlation terms) 

-I-(higher order terms)]. (9) 


In the above equation, the cross correlation terms and the 
higher order terms come from the fact that the brightness 
temperature is expressed as Eq. and they should appear 
even if the statistics of 5m , Sx-g^ and 5r, are completely Gaus¬ 
sian. In this sense, the first three terms in the above expres¬ 
sion, which are the auto-bispectra of 5m, and 5,,, should 
be corresponding to the intrinsic non-Gaussian features of 
these components and we focus on these auto-bispectra be¬ 
low. 

In Fig. O we plot the brightness temperature bispec¬ 
trum and the above auto-bispectra terms for equilateral type 
as functions of redshift (upper panels). From this figure, we 
can see that the total bispectra are mostly contributed from 


the auto-bispectra of the matter density field, the fluctua¬ 
tions of the spin temperature and the neutral fraction, which 
are expressed as the first three terms in Eq. m, for all red- 
shifts. For comparison, the evolution of power spectra is also 
shown (lower panels) and we find that the behavior of each 
component is very similar between bispectrum and power 
spectrum. Such a correspondence is highly non trivial, since 
the bispectrum and power spectrum reflect different aspects 
of the statistical properties of the fluctuations as we have 
mentioned. 

Let us try to interpret the behavior of bispectra, com¬ 
paring that of power spectra which w as detailed in our pre¬ 
vious work (IShimabukuro et al. 201^ . First, fluctuations in 
neutral hydrogen fraction appear when reionization begins 
and become dominant as reionization proceeds {z < 12). 
The dip at 2 ~ 14 corresponds to the redshift when the av¬ 
erage spin temperature becomes equal to the GMB temper¬ 
ature and the average brightness temperature 5Tb vanishes. 
This dip appears in the contribution of matter fluctuations 
for the same reason. Thus, this dip is independent of the 
properties of fluctuations and this is why both the power 
spectra and bispectra from and 5m have a dip at the 
same redshift. Note that we can also see the dip at z ~ 8 in 
the 21cm bispectrum at fc = O.lMpc”^. It would be due to 
sample variance coming from calculation of ionized bubbles 
because evolutions of ionized bubbles are different among 
other realizations. 

On the other hand, the spin temperature fluctuations 
are negligible at low redshifts {z < 10), because spin tem¬ 
perature is much higher than GMB temperature everywhere, 
that is, 7/ = 1 — T-y/Ts is very close to unity and indepen¬ 
dent of Ts- However they substantially contribute at higher 
redshifts {z > 14) and have two peaks at 2 ~ 15 and 2 25 

at large scales while the higher-redshift peak is much less 
noticeable at small scales. The dip at 2 ~ 23 is induced by 
the onset of X-ray heating. At higher redshifts (2 > 23), the 
spin temperature in dense region is lower than the average 
due to the WF effect, which couples the spin temperature to 
the kinetic temperature which is much lower than the GMB 
temperature. As a consequence, the probability distribution 
function (pdf) of spin temperature is negatively skewed at 
this epoch. Contrastingly, at lower redshifts (2 < 23), X- 
ray heating becomes effective for our parameter set and the 
spin temperature in dense region rises rapidly so that the 
skewness of the pdf changes its sign. Thus, at the onset of 
X-ray heating (2 ~ 23), the pdf becomes close to Gaus¬ 
sian with relatively small width, that is, the the skewness 
vanishes and the varianc e has a local minimum. Thus, in 
dShimabukuro et al. 201^ 1. it was suggested that the sign of 
skewness can be an indicator of the effectiveness of X-ray 
heating. The skewness and the variance have minima that 
coincide in redshift, similarly. It Is natural that the dips in 
the bispectrum and the power spectrum are coincident with 
each other. 

Before EoR, the bispectrum, as well as the power spec¬ 
trum, is mostly dominated by spin temperature fluctuations 
at large scales and will be a good probe of astrophysical 
effects such as the WF effect and X-ray heating. On the 
other hand, at small scales, the matter fluctuations are dom¬ 
inant and the bispectrum is of cosmological interest because 
it is induced by gravitational nonlinearity and, possibly. 
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z z 

Figure 5. Components of 21cm bispectrum for equilateral type: the brightness temperature (red), the contribution from eta (green), 
matter fluctuations (blue) and neutral hydrogen fraction (magenda). 


primordial non-Gaussianity (see, e.g., JScoccimarro 200d: 
iBernardeau et al. 20021 : Ijeong fc Komatsu 2009l 'l'). 

Next, we focus on the shape dependence of the total 
bispectrum and its components. Fixing fcs = 1.0 Mpc~^, we 
plot contours of the bispectra in (fci/fc 3 )-(fc 2 /fc 3 ) plain in 
top of Fig. [S] Note that we do not use the normalised bis¬ 
pectrum, kik^k^Blki, k 2 , ks), but the unnormalized bispec¬ 
trum, i3(fci, ^ 2 , fcs), here. We can see in what configuration 
of triangle the bispectra are strong. Here it should be noted 
that the triangle condition is not satisfied in the blank re¬ 
gion and that the contours are symmetric with respect to a 
line fci/fcs = fc2/fc3. 

At 2 = 10.05 when EoR has proceeded to some 
extent(a;i = 0.77), the total bispectrum is strong at folded 
and squeezed types. The contribution from neutral hydro¬ 
gen fraction fluctuations is dominant at these configurations, 
while matter fluctuation is dominant at equilateral type. At 
2 = 14.47, the dominant contribution comes from the spin 
temperature fluctuations and it is largest at squeezed type. 
The situation is similar at 2 = 27.03. At 2 = 20.23, both 
squeezed and folded type of the total bispectrum are strong. 
The contributions from matter and spin temperature fluc¬ 
tuations are comparable at these configurations, while the 
former is dominant at equilateral type. 

We also show the contour for fcs = 0.4 Mpc“^ in bot¬ 


tom of Fig|6] Compared with the case of ks = 1.0 Mpc“^, 
contributions from both matter and spin temperature fluc¬ 
tuations are significant at 2 = 20.23. On the other hand, 
we find that the contribution from fluctuations of neutral 
hydrogen fraction at 2 = 10.05 is clear compared with the 
case of ks = 1.0 Mpc“^. This helps us to subtract the in¬ 
formation on neutral hydrogen and it is better to see larger 
scales if we would like to know the information on neutral 
hydrogen fluctuations. 

4 DISCUSSION & SUMMARY 

In this paper, we investigated the 21cm bispectrum as a 
method to measure non-Gaussianity of brightness tempera¬ 
ture field. 

First, we have shown the scale-dependence of the 21cm 
bispectrum for the “equilateral-shape” at some redshifts. We 
found that the normalised 21cm bispectrum seems not to 
have any characteristic scale in 0.1 < fc/Mpc“^ < 1.0 for 
each redshift. For 2 = 10,15, 27, the normalised bispectrum 
is almost scale-invariant, while for 2 = 20 it has a scale- 
dependence as k^B oc k^. 

We have also shown the redshift evolution of the 21cm 
bispectrum with fixed k for three types of the shape in k- 
space. We found that the redshift evolution of the 21cm 
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Figure 6. (Top) Contours of the total bispectrum and its components in k\ jkz-ki jkz plane with fca = 1.0 Mpc (Bottom) Contours 
of the total bispectrum and its components in kxjk^-kijk^ plane with ^3 = 0.4 Mpc“^. 


bispectrum for the equilateral and folded shapes basically 
traces that of the 21cm power spectrum, but in case with 
the squeezed shape, we could see a different behavior and it 
can be understood by considering the coupling between the 
large- and small-scale modes. 

Then, we studied the 21cm bispectra by decomposing it 
into the contributions from the matter density field, the fluc¬ 
tuations in the spin temperature and the neutral fraction. 
From the redshift evolution, we found the dominant com¬ 
ponent at each redshift and scale. We also show the shape 
dependence of each component and compared it with that 
of total 21cm bispectrum. The shape dependence of each 
component looks similar to each other, but a slight differ¬ 


ence also exists. Hence, by future precise observation it is 
expected that we would obtain the information about the 
non-Gaussian feature of these components separately. 


As far as the matter bispectrum is concerned, there 
have been a lot of works which discuss the shape-dependence 
by using the second order perturbation theory and also 
numerical N-body simulation. Although most of the works 
focus on the matter bispectrum at the lower redshift 
[z ~ 1.0) or higher redshift (30 < 2 5, 100 in the dark 
age) _ (e.g., Js coccimarro 2 OOOI: iBernardeau et al. 200: 


Jeong fc Komatsu 20091 : ISefusatti et al. 2010l : [1|; 

Munoz et al. 20l5h ). we can find that the shape of the 
matter bispectrum in our result is basically consistent with 
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these previous works by extrapolation. By using the second 
order perturbation theory, the matter bispectrum can be 
expressed as Bm{ki,k 2 ,k-i) oc Pm{ki)Pm(k 2 ) + 2 perms, 
with Pm{k) being the matter power spectrum. For 
0.1 Mpc“^ < k, the matter power spectrum behaves 
as oc Based on this fact and the isosceles ansatz 

(fci = ^2 = fc = aks), we have Bm{k, a) oc + 
and hence the unnormalised matter bispectrum becomes 
larger as a increases is the largest in the squeezed shape 
lljeong fc Komatsu 2009l i. 

Based on the above discussion about the matter bis¬ 
pectrum, let us revisit the behavior of the 21cm bispectrum 
as shown in Fig. [T] As we have mentioned, in contrast to 
other redshifts, for z = 20 the normalised 21cm bispectrum 
has a scale-dependence as k^ in the equilateral shape. From 
Fig. [5] we find that at 2 = 20 the matter contribution rela¬ 
tively dominates over the 21cm bispectrum, and hence the 
behavior of the 21cm bispectrum is expected to trace that of 
the matter bispectrum at this redshift. Based on the expres¬ 
sion obtained from the second order perturbation theory, 
the scale-dependence of the normalised matter bispectrum 
can be estimated as oc fc® x Hence, the be¬ 

havior of the matter bispectrum could explain that of the 
21cm bispectrum at 2 = 20. 

Note that the Zel’dovich approximation is used to 
solve density evolution in 21cmFAST. This approximation 
is imperfect to calculate the bispectrum and it is desir¬ 
able to use N-body simulations to estimate the bispec¬ 
trum accurately. However, at scales less than k ^ 5Mpc~^ 
at 2 = 7 — 20 , density evolution with the Zel’dovich 
approximat ion coincides with t hat obtained by N-body 
simulations dMesinger et al. 201lli . 

Naively, the spatial distributions of the spin tempera¬ 
ture and neutral fraction should be considered as tracers of 
matter density field, that is, they could be treated equally 
with those of the halos and galaxies. Based on this consider¬ 
ation, the simplest way to express the spatial distributions 
of the spin temperature and neutral fraction is introducing 
a bias parameter, such as Si oc biSm {i = rj and xh). If 
such a bias parameter is scale-independent, the behaviors of 
the bispectra of the spin temperature and neutral fraction 
are completely the same as that of the matter density held. 
However, we can see slight differences between these com¬ 
ponents in Fig. |6] and also[T] Hence, we expect that the bias 
parameter should have non-trivial scale-dependence due to 
the non-linear or non-local transfer from the matter density 
held to the spin temperature and neutral fraction. We need 
to investigate this issue more deeply in future work. 

Here, we briehy discuss the cross correlation higher or¬ 
der terms in eq.(|9]). Although we only plot auto correlation 
term in Fig[S] neglected cross correlation and higher order 
terms are of the same order as the auto correlation terms. 
We show the ratio between the sum of three component 
terms and total 21cm bispectrum as function of redshift in 
Figl3 As you can see from this hgure, this ratio is not unity 
since we do not include the cross correlation and higher or¬ 
der terms, and hence, in practice, these terms can not be 
neglected. However, in our study, we focus on only auto 
correlation terms in order to study contribution from each 
component. Although we expect that what cross correlation 
terms are effective at each epoch from Fig|6] (e.g. cross cor¬ 
relation between fluctuations neutral hydrogen fraction and 


matter density might become effective at the EoR (2 ~ 10)), 
detailed study is our future work. 

The detectabili ty of bispectrum is of critical interest. In 
our previous work llYoshiura et al. 2015^ . we estimated the 
signal-to-noise ratio of bispectrum, developing a formalism 
to calculate the bispectrum contributed from thermal noise. 
We find that the SKAl has enough sensitivity for both equi¬ 
lateral and isosceles(A' = |fci| = \k 2 \,k = Ifcsl with K=0.06 
Mpc“^) triangles for k < 0.3 at 2=8-17, while LOFAR will 
have sensitivity for the peaks of the bispectrum as a function 
of redshift. Actually, galactic and extragalactic foreground 
will be a serious obstacle just as in the case of power spec¬ 
trum and should be studied in detail. 

There are some other approaches to measure non- 
Gaussianities in the brightness temperature field. For ex¬ 
ample, some studies focus on topological structure of 
brightness temperature field such as Minkowski fu nctionals 
llGleser et al.2()0^ : iLee et al. 20081 : iHong et al. 20141 : [3) . This 
method is complementary to higher order statistics. 


APPENDIX A 


In this section, we derive the relation between bispectrum 
and skewness. When we define the fluctuation of brightness 
temperature STi, as 5 = {STi, — 5Tb)/STb, the skewness of 
brightness temperature is defined by 

■I ^ _ 

^ = j;^J2(STb,i-STbf (10) 

i = l 

^ {STb)^ ^3 

N ^ 

= mf{s^). 


Here N is the number of pixel. The definition of three- 
point correlation function for brightness temperature ^ is 
expressed by 5(ri, r 2 ) = ((5(x)5(x-|-ri)(5(x-f r 2 )). We can 
connect this three-point correlation function with ensemble 
average of <5® such as (5®) = ^(0,0) . We know the relation 
between correlation function and bispectrum from Wiener- 
Khintchine relation described by 


C(ri,r2) 



d k2 qki.ri+k2-r2) 
(27r)3® 


i3(ki, k2, 


-ki-k2). 

( 11 ) 


Therefore, ^(0,0) can be expressed with bispectrum by 


= / / ^S(k„k2,-kr-k2). (12) 

By using eauation jllll and (I12II . we can derive the relation 
between skewness and bispectrum as. 




(27! 


rB(ki, k2,-ki — k2). (13) 
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Figure 7. This shows the fractional difference between sum of three component terms and total 21cm bispectrum at A: = O.lMpc ^{left) 
and k = 1.0Mpc~^ (right). Here, we plot the (Bg^ + + B^h)/B gTj^ 
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